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We use analytical arguments and large-scale Monte Carlo calculations to investigate the nature of the phase 
transitions between distinct complex superfluid phases in a two-component Bose-Einstein condensate when 
a non-dissipative drag between the two components is being varied. We focus on understanding the role of 
topological defects in various phase transitions and develop vortex-matter arguments allowing an analytical 
description of the phase diagram. We find the behavior of fluctuation induced vortex matter to be much more 
complex and substantially different from that of single-component superfluids. We propose and investigate 
numerically a novel drag-induced "preemptive vortex loop proliferation" transition. Such a transition may be a 
quite generic feature in many multicomponent systems where symmetry is restored by a gas of several kinds of 
competing vortex loops. 
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I. INTRODUCTION 

Natural generalizations of many superfluid phenomena are 
possible in mixtures of independently conserved multicompo- 
nent Bose-Einstein condensates with intercomponent cuiTent- 
current interactions. The topic was first investigated in the 
context of *He — ^He mixturesi^, where it is possible to at- 
tain only a limited range of parameters. The recent progress in 
atomic Bose-Einstein condensates (BEC) has made it possi- 
ble to access a much wider range of regimes and explore novel 
superfluid phases which can arise in such mixtures. For this 
reason, there has been much interest in a generic example of 
an interacting BEC mixture, namely a C/(l) x t/(l)-symmetric 
system with cuiTent-cuiTent interactions. One of the novel as- 
pects of the superfluid physics in such a system is the pos- 
sibility of a phase transition at a sufficiently strong current- 
current interaction to a state of paired superfluid where the 
only broken U{1) symmetry is associated with order only in 
the phase sum^. The other discussed example (which does 
not fall within the framework of Galilean-invariance based 
argument^) is a phase transition for bosons on an optical lat- 
tice to a state where one species of bosons pair with holes of 
the other species, and thereby retaining order only in the phase 
differenced^. 

These transitions were investigated numerically in great 
detail in the J-current model in Ref. 6 using the worm- 
algorithm^. This numerical study, combined with mean-field 
arguments, revealed the interesting fact that with increasing 
cuiTent-cuiTent interaction, the usual second-order superfluid 
phase transition is altered to a first order phase transition^. 
In the free energy functional, the current-current interaction 
is consistent with U{1) x U{1) symmetry and the transition 
should therefore be associated with a proliferation of interact- 
ing vortex loops where all vortex-loop segments of the system 
interact with each other through a Coulomb potential. Ex- 
isting theories of proliferation of such defects, however, al- 
ways lead to a second-order superfluid phase transition**. This 
indicates that in this system we are faced with a novel sce- 
nario for thermally driven spontaneous vortex-loop prolifer- 



ation, the detailed investigation of which is the goal of the 
present work. 

To describe the behavior of the system undergoing these 
phase transitions as proliferation of vortex loops in a two- 
component condensate, we propose a scenario of a "pre- 
emptive vortex-loop proliferation". This scenario in par- 
ticular allows us to estimate the characteristic critical cou- 
plings (or equivalently, critical temperatures) and provides a 
vortex-matter based picture of the transitions in the most in- 
teresting part of the phase diagram, from a state with broken 
U{l)xU{l) symmetry into a paired superfluid state and a sub- 
sequent transition into a normal state. To find numerical back- 
ing for the preemptive vortex-loop proliferation scenario, we 
perform a large-scale Monte Carlo (MC) calculation of vortex 
matter in the interacting BEC mixture using a representation 
in terms of the phases of the ordering fields of the conden- 
sates. This numerical approach allows us to study directly vor- 
tex matter and therefore may be viewed as complementary to 
the worm-algorithm based approach in Refs. .3.6>,12t The in- 
sight which we obtain from Monte Carlo calculations on vor- 
tex matter may also shed light on how the Andreev-Bashkin 
effect' modifies the vortex-matter phase transition predicted 
for the liquid metallic state of hydrogen^ 

Finally, we remark that the problem of multicomponent 
vortex-loop proliferation has a quite generic character, since it 
is also related to a wide spectrum of phase transitions in other 
systems. An example is represented by individually conserved 
electrically charged condensates that communicate with each 
other only via a fluctuating gauge fie\(fd2ilLiZil2JdiH, More- 
over, a related problem arises in three-dimensional generaliza- 
tions of phase transitions discussed recently for certain planar 
spin-1 condensate a'^''^ . 



n. THE MODEL 

We consider a mixture of Bose-Einstein condensates with 
U{1) X U{1) symmetry and current-current interaction. This 
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system in the hydrodynamic limit is described byi 

^ ^\ J {{pi - pd)vj + {p2 - pd)vl + 2pdVi ■ V2} 
= ^ J^dr{pivl+ p2vl~ pd{vi-V2f} , (1) 

where Vi — hV9i/mi. The last term describes a current- 
current interaction' (for its detailed microscopic derivation, 
see Ref. ITsI) . The microscopic origin of the non-dissipative 
drag can for example be elastic inter-component scattering 
due to van der Waals forces between the charge-neutral atoms 
in the system.'^, or can also originate from a lattice-'-^. This 
coupling is consistent with U{1) x U{1) symmetry and thus 
is very different from the symmetry breaking intercompo- 
nent Josephson-coupling, which is a singular perturbation. A 
drag term is perturbatively irrelevant and a critical strength is 
needed to change the zero-drag physics of the problem, due 
to the extra two gradients in the coupling between the two 
phases. 

The discrete model as such may also have a physical real- 
ization in terms of a Bose-Einstein condensate on an optical 
lattice^. In the latter case, a particularly wide range of both 
positive and negative pd can be accessed^. The parameter pd 
is a superfluid density of one condensate carried by the su- 
perfluid velocity of the other as follows from the equations of 
motion'. 



Jl = (Pl - Pd)vi + PdV2, 
32 = (P2 - Pd)v2 + PdVl- 



(2) 

(3) 



Symmetry-restoring phase transitions in this system are asso- 
ciated with proliferation of thermally excited topological de- 
fects, namely vortex loops**. In what follows, we denote vor- 
tices in the two-component condensate by a pair of integers 
corresponding to the winding of the phases in each of the con- 
densates 



(A6I1 = 27rni, A6'2 = 27rn2) = (^1,^2). 



(4) 



The current-current interaction 2pdVi ■ V2 introduces a bias 
for counter-directed currents when pd is positive . Indeed, this 
term introduces an attractive Coulomb interaction between 
(±1,0) and (0, =f1) vortices. The coefficients pi,p2 and pd 
must satisfy the relation 



Pd < 



P1P2 
Pi+ P2' 



(5) 



for stability. This puts an absolute upper bound on the amount 
of drag in the system that can be considered physical. In the 
phase diagrams to be presented below, we denote as gray (for- 
bidden) those areas which cover the sets of parameters that 
violate the above inequality. 



III. ENERGY SCALES ASSOCIATED WITH BARE 
STIFFNESSES 

Let us begin by a straightforward examination of the en- 
ergy scales of the problem. In what follows, we set the masses 



of the condensates equal and absorb h/m in the definition of 
p. We focus on the pd > case. In what follows, we de- 
note expressions for phase stiffnesses for various topological 
excitations as J(i.j), (the index refers to corresponding 
topological defect). The explicit expressions are given by 



•^(1,0) = Pi - Pd, 



J{0,1) = P2 



JilS) = Pl+ P2, 



J(l^_l) = Pi + P2 - 4pd. 



(6) 



(7) 



(8) 



(9) 



Let us denote the critical stiffness of the 3DXF- Villain 
model''' by 



Pc « 0.33 



(10) 



Then if we neglect any interactions between different species 
of vortices, naive estimates of the lines where various vortex 
modes would proliferate, are given as follows. 

(1, 0)-vortices proliferate from an ordered background 
along a line defined by pi — pd — Pc- 

(0, l)-vortices proliferate from an ordered background 
along the lines defined by p2 — Pd = Pc- 

(1,-1) vortices would proliferate from an ordered back- 
ground along a line defined by pi + p2 — ipd — Pc- 

Proliferation of (1, 1) is irrelevant because of the above 
types of topological excitations always proliferate (and thus 
restore symmetry) before (1, 1) vortices, when pd > 0. 

Below we show that this naive energy-scale based picture 
is not correct. 



IV. PHASE DIAGRAM, EQUAL STIFFNESSES 

The simplest case is where the bare phase stiffnesses of 
each component is equal, so we begin by considering that case 
first. 



A. Continuous phase transitions in limiting cases 

The character of the vortex-loop proliferation transition can 
readily be understood in two limiting cases, by mapping the 
system to a single component model yielding standard second 
order phase transitions. 

One limit is the trivial limit pd 0, when the system is de- 
scribed by two independent XY models undergoing a second 
order phase transition from U{1) x J7(l)toa symmetric state. 
Indeed, in this limit there is no energetic or entropic advantage 
in restoring order by composite topological defects. 

Another limit which is fairly simple to understand, follows 
from the fact that by increasing pd, the stiffness of (1,-1)- 
composite defects can be made arbitrarily much smaller than 
the stiffnesses for (1, 0) and (0, 1) defects. This is the limit 
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where 2p — Apd ~ Pc < {p — Pd) and thus the vortex loop 
(1, —1) costs Uttle energy to excite, while (1, 0) and (0, 1) ef- 
fectively are frozen out. Physically, this also means that in this 
limit it is energetically costly to split a composite (1,-1) de- 
fect into a pair of individual vortices, and therefore one may 
neglect its composite nature and map the system onto a 1- 
component 3DXY model undergoing a phase transition at 
J(i _i) ~ 2p ~ Apd = Pc- Because (1,-1) vortices can- 
not disorder the phase sum, this continuous phase transition 
is associated with going from a U{1) x U{1) state to a state 
with U{1) symmetry associated with order in the phase sum, 
which is the "paired superfluid phase" in Ref. 3. 

Let us now consider the other regimes which occur in the 
case pi = p2 = P case. For some regimes another representa- 
tion of Eq. ([T]l will be useful, namely 

+ ^[V(0i+02)]'}. (11) 

This form of the energy is particularly useful when we want 
to discuss the vortex matter of the remaining superfluid com- 
ponent in the background of proliferated composite vortices. 
We next proceed to discussing this situation. 

B. Phase transitions in a nontrivial vortex gas background 

A deviation from the vortex proliferation based on the naive 
energy scales scenario is manifested in the transition to a fully 
symmetric state in the regime J{i.^i) < J{i,o) = ■^(o.i)' i-^- 
2/5 — Apd < Pc < p — Pd- To understand how this takes 
place, we should understand how the background of prolifer- 
ated (1,-1) vortices affects (1,0) and (0,1) vortices. This 
can be explained from the separation of variables in Eq. ( fTTT ). 
The spontaneous proliferation of (1,-1) vortices leaves the 
remaining broken symmetry only in the second term. The 
corresponding remaining phase stiffness is that of a "clapping 
mode" associated with a response to varying the phase sum. 
The stiffness of the clapping mode is destroyed by prolifer- 
ation of the cheapest topological defects with a winding in 
the phase sum. These defects are individual vortices (1,0) or 
(0,1). The separation of variables Eq. ( fTTl i suggests that the 
background of proliferated (1,-1) vortices destroys the phase 
stiffness in the first term and thus only the second term deter- 
mines the effective stiffness of (1, 0) or (0, 1) vortices. Their 
stiffness is therefore reduced compared to the bare stiffness in 
the naive energy scale argument. The new effective stiffness 
is >^(i,o) — J{i.o) = and thus it suggests that the system 
undergoes a phase transition to a fully symmetric state at 

1 

2P = Pc- (12) 

Note that from this argument, it follows that the proliferation 
of (1,0) or (0, 1) vortices in the background of proliferated 
(1,-1) vortices is determined by p only. This is testable in 
MC calculations, and we report on it below. 



C. Preemptive phase transition 



Now consider the most interesting regime where the line 
defined by the relation J(i,o) = '^(o,i) — P ~ Pd = Pc inter- 
sects the line defined by the relation J(i_i) = 2p~Apd — Pc- 
We denote the intersection point derived from the naive energy 
scale-based argument by {pj,pdi) — {Spc/^, pc/2). Con- 
sider the regime slightly above the point pi — 3pc/2 (i.e. 
p — Pi + S and pd ~ Pdi + 6/2). Then, from Eq. (fTTI ) 
we conclude that although the phase transition is indeed ini- 
tiated by proliferation of the lowest-in-energy topological de- 
fects ((1, —1) in this regime), the remaining stiffness for (1,0) 
and (0, 1) excitations p/2 w pi/2 — 3pc/4 (which can be 
read off from the second term in Eq. (fTTT)). is actually less 
than Pc- Hence, the vortices (1, 0) and (0, 1) cannot remain 
confined once (1,-1) are proliferated. Therefore, from the 
separation of variables we may draw the conclusion that the 
simple energy-scale based picture underestimates the critical 
stiffnesses. More importantly, away from the limiting cases, 
the process is cooperative and hence proliferation of compos- 
ite defects may trigger proliferation of individual vortices at 
a critical stiffness where arguments based on energy scales 
alone would predict that the individual vortex loops remain 
confined. Thus, with respect to (1,0) and (0,1) vortices, 
we are dealing with a "preemptive" vortex-loop proliferation 
scenario, triggered by the interaction with vortices in a dif- 
ferent sector of the model. In the case where the energy of 
(1, —1) vortices is almost the same as that of (1, 0) and (0, 1) 
vortices there is only one transition where by the same ar- 
guments both types of topological defects assist each other 
in restoring symmetry via a single phase transition. Numer- 
ical calculations which we report in the second part of this 
paper confirm this behavior of vortex matter Importantly, 
whenever we observed this behavior, the phase transition was 
first order within the resolution limits of our MC calculations. 
The region of the phase diagram showing first order transi- 
tions in our computations, appears to be consistent with the 
findings in the J-current model^ with the same symmetry, 
though in our case the microscopic physics is different. Note 
that this scenario is substantially different from the continu- 
ous loop-proliferation transition invariably encountered in a 
single-component model^ii^. 

Fig. [T] summarizes the new estimates for the lines of vor- 
tex proliferation which follow from the separation of variable 
argument Eq. ( fTTT i. They are given by three different regimes. 

(1, 0)- and (0, l)-vortices proliferate from an ordered back- 
ground along a line defined hy p — pd = Pc (solid red line in 
Fig.D, 

(1,-1) vortices proliferate from an ordered background 
along a line defined by 2p — Apd — Pc (dashed blue fine in 
Fig.[D, 

(1,0)- and (0, l)-vortices proliferate from a background 
of proliferated (1. — l)-vortices at p/2 — pc (dashed-dotted 
black line in Fig. [T]). 
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FIG. 1 : (Color online) Schematic phase diagram which follows from 
separations of variables for the model in Eqs. (TJ and dllb . pi — 
p2 = p. The gray-shaded area is the forbidden regime pd > p/2. 
The lines separating the various regions are: a) Solid (red) line: Pro- 
liferation of (1, 0) and (0, l)-vortices from an ordered background 
[i.e. large p), along the line pd = p — Pc- b) Dashed (blue) line: Pro- 
liferation of (1, — l)-vortices from an ordered background, along the 
line Pd = p/2 — pc/4. c) Dashed-dotted vertical (black) line: Pro- 
liferation of (1, 0) and (0, l)-vortices from a background of prolif- 
erated (1, — l)-vortices, along the line p = 2pc- This vortex-matter 
phase diagram has the same topology as that obtained from the J- 
current model^^. I: (7(1) x (7(1); 11: L''(l)-symmetry in the phase 
sum; and III: a fully symmetric case. 



V. PHASE DIAGRAM, UNEQUAL STIFFNESSES 

We next generalize the above qualitative considerations to 
the case of unequal stiffnesses pi ^ P2- We will use the nota- 
tion that pi = api : pi = p. So that the coefficient a is a mea- 
sure of the disparity of the stiffnesses. Since the non-triviality 
of the phase diagram of this model is associated with the pos- 
sibility of tuning the energy of composite (1,-1) defects to be 
less than (or comparable to) the energy of individual vortices 
(1, 0) and (0, 1), we use the same strategy as in the previous 
section to analyze this model. That is, we need to separate the 
variables by extracting all the stiffness terms which are unaf- 
fected by proliferation of (1,-1) vortices. The corresponding 
part of the free energy functional therefore should depend on 
gradients of the phase sum only. Separating the variables in 
such a way we arrive at the following representation of the 
model 



F=- I dr 



ap^ - (1 + a)ppd 



[V (01+02)]' 



(1 + a)p - Apd 



{l + a)p-Apd 
[(p-2prf)V0i -iap-2pd)V02? 



(13) 



Notice the asymmetric phase weights in the second term in 
contrast to the symmetric separation of variables in Section 
HV] The asymmetry of the problem is also seen if we consider 
negative pd which would result in decreasing the energy of 
(1, 1) vortices compared to (1, 0) and (0, 1) vortices. Negative 
Pd may be easily realized in Bose-Einstein condensates on an 
optical lattice, and we consider this possibility in Appendix B. 



In the following qualitative discussion in this section, we 
consider only a positive pd and without loss of generality, we 
assume that a > 1. We start by going through the same en- 
ergetics as we did for the case pi ~ P2- First of all, the con- 
dition for stability Eq. (|5]l, now reads pd < a.p/{l + a). 
The proliferation of (0, 1) -vortices from an ordered back- 
ground is now determined by the condition ap — pd — Pc 
or equivalently pd ~ ap — Pc while that of (1, 0)-vortices 
is determined by the condition p — pd — Pc, or equivalently 
Pd = p — Pc- These lines now differ from each other, in con- 
trast to the case pi — p2, and hence there will be one ad- 
ditional region in the phase diagram. This follows, since at 
Pd — 0, the phase transitions in the model are expected to 
be two non-degenerate vortex-loop proliferation transitions 
in the 3DXy-universality class, with a regime with ordering 
only in one phase separating them. Moreover, the spontaneous 
proliferation of (1, — l)-vortices from an ordered background 
is now determined by the condition (to be read off from the 
second term in Eq. ( fT3] l) {1 + a)p — Apd = Pc or equivalently 
Pd = [(1 + a)p — Pc]/4:. These expressions reduce to those 
that were discussed in Section lTV Cl for the case a = 1. 

Next, we proceed to investigate the condition for prolifer- 
ation of (0, 1)- or (1, 0)— vortices in a background of prolif- 
erated (1, — l)-vortices. We thus assume (an assumption that 
will be checked numerically in the second part of the paper) 
that the coefficient of the stiffness associated with the second 
term in Eq. ( fT3] ) has renormalized to zero. Then, the first term 
accounts for the only phase stiffness remaining in the system. 
Thus, the effective model becomes 



^(1,-1) 



= ^/dr^^f^-(i±4^[V(0i+02 
2 (l + ajp-Apd 



(14) 



Note the rather surprising fact that, provided the compos- 
ite vortices (1,-1) have proliferated, the (1,0)- and (0, 1)- 
vortices enter the effective model on equal grounds even if the 
bare phase stiffnesses for these differ The origin of this fact is 
that (1,-1) composite defects have an asymmetric effect on 
the partial reduction of the bare phase stiffnesses of the indi- 
vidual vortices. Thus (1, 0)- or (0, 1) vortices will participate 
on equal grounds in the restoration of the remaining symme- 
try. Based on the above conjectures we obtain the condition 
for proliferation of (1, 0)- or (0, 1) vortices in the background 
of proliferated (1,-1) loops 



gp^ - (1 + a)ppd 
{l + a)p-Apd 



(15) 



Observe that in contrast to the similar condition Eq. (fTSl i for 
the case of equal stiffnesses in Section lTVCl when a 1, pd 
no longer drops out of this relation. The explicit relation is 



Pd 



ap^ - (1 + a)ppc 
(1 + a)p - Apc 



a ^ 1, 



(16) 



which is seen to approach pd — > ap/(l + a) from below 
as p becomes large, i.e. the proliferation line approaches the 
forbidden parameter region from below as p oo. The 
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line of proliferation under discussion, namely the proliferation 
of (1,0)- and (0, l)-vortices in the background of prolifer- 
ated (1, — l)-vortices, only comes into play above the dashed 
(blue) line separating phases I and II in Fig. |2] This is when 
the composite vortices are actually proliferated. We therefore 
only plot the line in this regime, and this is the dashed-dotted 
(black) line given in Fig. |2] 

In Fig. 12] the solid (red) lines are the lines of proliferation 
of (1, 0)- and (0, l)-vortices from an ordered background. At 
Pd = 0, they emanate linearly from p — Pc/ct for (1,0)- 
vortices growing as ap, and from p = pc for (0, l)-vortices 
growing as p. The dashed (blue) line represents the line of 
proliferation of (1, — l)-vortices from an ordered background. 
It emanates at p^; = from p = pc/{l + a), growing as 
[(1 + a)p]/4:. The dashed-dotted (black) line represents the 
line across which the effective stiffness of the clapping mode 
01 + 02 vanishes through the proliferation of individual vor- 
tices (1, 0) or (0, 1). The lines are seen to divide the phase 
diagram into four distinct regions, namely I) the completely 
ordered state, II) the partially ordered state with proliferated 
(1, — l)-vortices and confined individual (1,0)- and (0,1)- 
vortices. III) the completely disordered state with proliferated 
individual vortices, and IV) the partially ordered state with 
confined (0, l)-vortices and proliferated (1, 0)-vortices. Re- 
gions II and IV are therefore two distinct partially ordered 
states with one broken U (l)-symmetry in each case. 



Pd/Pc 




FIG. 2: (Color online) Schematic phase diagram with regions I,II, III, 
and IV for the model in Eqs. l[T) and l ll3t with p2 — api > pi — p. 
For the purposes of illustration, we have taken a = 1.2. The gray 
shaded are is the forbidden parameter regime pd > ap/{l + a). 
The lines separating the various regions are obtained as follows: 
a) Solid (red) lines: Proliferation of (0, l)-vortices from an ordered 
background, along the line pd = ap — pc, as well as prolifera- 
tion of (0, l)-vortices from an ordered background, along the line 
Pd ~ p — Pc- b) Dashed (blue) /;«e: Proliferation of (1, —l)-vortices 
from an ordered background, along the line Pd — (l + a)p/4 — pc/4. 
c) Dashed-dotted (black) line: Line of proliferation of individual vor- 
tices (1, 0) or (0, 1) in a background of proliferated vortices (1, —1), 
given by Eq. [16] When we cross this line from right to left, passing 
from region II to III, the stiffness of the clapping mode 6i + 62 is 
destroyed by the proliferation of individual vortices. 



A. Preemptive scenario 

We next discuss the preemptive scenario for vortex-loop 
proliferation for the more general case pi ^ p2, largely fol- 
lowing the line of reasoning in Section IIVCI It turns out 
that the physics is quite rich and markedly different from 
the single-component case, which is rather surprising given 
the simplicity of the coupling term between the two conden- 
sates, cf. Eq. ([T]). Hence, consider the intersection point 
where the line of proliferation of (1, — l)-vortices from an or- 
dered background is intersected by the line of proliferation 
of (1, 0)-vortices from an ordered background. This intersec- 
tion takes place at {pj,pdi) = (3pc/ (3 — a), apc/{3 — a)). 
Consider now a point slightly above the intersection point 
above the line defined by the relation {pi + 6, pdi + 5d) = 
{3pc/ (3 — a) + S, apc/ (3 — a) + (1 + a)6/4), where compos- 
ite vortices are proliferated. The remaining stiffness for the 
clapping mode, pciap, is given by Eq. (O 



Pclap 



gp^ - (1 + a)ppd 
{l + a)p-Apd 



(17) 



The question is now whether proliferation of (1, —1) vortices 
can trigger a preemptive proliferation of (1, 0) and (0, 1) vor- 
tices. By evaluating PcUpiPi Pd^ct) at the intersection point 
{pi, Pdi) between proliferation of individual vortices in an or- 
dered background and proliferation of composite vortices in 
an ordered background, the issue is if Pc\ap{pi, Pdi, ct) < Pc, 
(i.e. if this estimate yields a situation that upon proliferation 
of composite vortices the individual vortices no longer have 
enough stiffness remaining to stay condensed). If this is the 
case, then our estimates will indicate a preemptive vortex-loop 
proliferation, following the same line of reasoning as was used 
in Section lTV CK to be checked in Monte Carlo calculations in 
the second part of the paper). 



VI. WEIGHTED PHASE SUM ORDER 

It has been observed in the past that in the drag problem Eq. 
(|T]i, the vortices of the type (1, —n) with n > 1 can become 
energetically cheapesti^. Let us apply the separation of vari- 
ables method to estimate analytically the position and drag 
dependence of the transition lines in the phase diagram when 
(1, ~n) -types of defects are relevant as well as to describe 
how vortex matter drives transitions from partially ordered to 
fully symmetric states in these cases. The accuracy of this 
method will be checked numerically in the second part of the 
paper 

Consider p2 < pi and pd > 0. First, one should ex- 
amine for which ratio of the bare stiffnesses p2 / pi does the 
system prefer to proliferate composite (1, — n — 1) vortices 
rather then (1, —n). The conditions when the energy for an 
(1, — ri — 1) excitation is less then that of an (1, — n) excitation 
can be found as follows. The phase stiffness associated with 
an (1, —n) excitation is J(i^_„) = pi + n^p2 ^ (1 + "■)^ Pd- 
Hence one finds that the inequality pi + {n + 1)^P2 — (n + 



2ypd < Pi 



n?P2 



[n + 1) Pd must be satisfied if the 
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system is to prefer proliferating (1, —n — 1) vortices in an 
ordered background instead of proliferating (1, — ri) vortices. 
Combined with the constraint pd < piP2/{pi + P2) on pd we 

which gives 



find |^P2 < Pd< 



P1P2 
P1+P2 



Pi n 



1 



1/2- 



(18) 



This condition is illustrated in Table |T] From this, it follows 
that for P2//O1 < 2/3, it is energetically less costly to excite 
(1. —2) vortices rather than (1, —1) vortices for sufficiently 
large value of pd- 



TABLE I: This table shows the condition for the ratio between the 
bare stiffnesses, when we assume that p2 < pi and pd > 0, for the 
system to proliferate a given composite vortex. 



Composite vortex 


Condition 


(1,-1) 


2/3 < P2/P1 < 1 


(1,-2) 


2/5 < P2/P1 < 2/3 


(1,-3) 


2/7 < pa/pi < 2/5 



For such regimes the proper separation of variables is 



1 



P1P2 - Pd {Pi + P2) . ™ , .2 



F = - dr 

2 .L \,pi+n'^p2-(l + nypd 



1 



pi + n'^p2 - (1 + nY 



(pi - (1 + n)pd) VOi 



P2-il + -)pd)y02 
n 



(19) 



were n is an integer. This separation of variables is performed 
in order to extract the part of the free energy which is unaf- 
fected by (1, —n) winding in the phases. Thus, upon prolifer- 
ation of (1, ~n) vortices the system enters a phase with order 
in the weighted phase sum nOi + O2 (while individual phases 
are disordered). The effective phase stiffness which will re- 
main in the system is given by 



(1,-n) 



- /dr- 



P1P2 - Pd jpi + P2) 
Pi +7i2p2 - (1 + nf Pd 



{nVOi + V6'2) (20) 



In contrast to the case considered in previous sections, here 
the individual phases do not participate on equal grounds af- 
ter proliferation of (1, — n) vortices because one of the phases 
has a factor n and is therefore more expensive to fluctuate. 
Nonetheless, there are several types of topological defects 
which can contribute on equal grounds to restore the remain- 
ing symmetry. In Fig. [3] we plot Pciapipi, Pdi,a)/pc as a 
function of a. 

For definiteness, we next consider in detail the case n ~ 2. 
Then, the cheapest topological defect with which to restore 
the symmetry in Eq. (|20] i is given a doublet of an elemen- 
tary vortex (0, 1) and a composite vortex (1, —1) which is of 
lower order than the vortex (1, —2) which drives the system 




a 



FIG. 3: (Color online) Plot of Pciap(p/, pdi,<^) I pc as a function of 
a. Note that the system is symmetric around 01 = \. The dashed 
(black) line in the above figure shows that the remaining stiffness 
of the clapping mode (in the background of proliferated compos- 
ite (1, — l)-vortices) is less than the critical coupling for vortex- 
loop proliferation in a parameter regime 2/3 < a < 3/2. This 
is the parameter regime where it is correct to limit oneself to the 
sector where the composite proliferated background vortices are of 
the type (1,-1). For a < 2/3, the composite proliferated back- 
ground vortices are of type (— n, 1), while for a > 3/2, the com- 
posite proliferated background vortices are of type (1, —n). The 
fact that Pcia.p{pi , Pdi,Q) < pc indicates that we are in a param- 
eter regime where the full restoration of (7(1) x I7(l)-symmetry 
proceeds from a preemptive vortex-loop proliferation phase transi- 
tion, as explained in the text. The solid (red) line in the above figure 
shows pciap (p/, Pdi, q) for a < 2/3 in a regime where (2, — 1)- 
vortices trigger preemptive proliferation of all topological defects 
and a > 3/2 where (1, — 2)-vortices initiate the phase transition 
into fully symmetric state. 



into the partially ordered state ( |20l i. The transition back to a 
fully symmetric state then takes place when 



P1P2 - Pd jpi + P2) 
Pi + n'^P2 - (1 + Pd 



Pc- 



(21) 



In our MC calculations, which we report below, we check this 
dependence. 

Before we proceed to the Monte Carlo calculations, we re- 
mark on the accuracy of the estimates of the location of the 
phase-transition lines based on separation of variables. The 
location of the phase-transition lines based on the above argu- 
ments, have corrections in the regimes of the phase diagram 
where several such lines split. This is because in the vicinity 
of such splitting points, the energy scales associated with vari- 
ous types of topological defects are not well separated. Hence, 
energetically next-to-cheapest excitations could participate in 
the depletion of the phase stiffness. The above arguments be- 
come more accurate as we move away from splitting points. 
However, they underestimate critical stiffnesses near splitting 
points. Below, we perform Monte Carlo simulations to study 
the least analytically tractable region near the splitting points. 
We find that even near the splitting points, the separation-of- 
variables based argument is quite accurate. 
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VII. MONTE CARLO CALCULATIONS 

We next proceed to presenting our numerical results based 
on large-scale Monte Carlo calculations, for which we need 
to define our continuum model on a numerical lattice. Al- 
ternatively, we may view it as a physical realization of a 2- 
component Bose-Einstein condensate on an optical lattice, 
as alluded to above. Providing a faithful lattice representa- 
tion of the continuum model Eq. ([T]i using phase variables 
is not straightforward, as some of the schemes for formulat- 



ing the theory on a lattice, which are standard in the single- 
component case, introduce subtle artifacts when the current- 
current interaction between two condensates is discretize. It 
turns out that a study of the vortex physics in a lattice repre- 
sentation of the model Eq. ([Til is best facilitated by the so- 
called Villain approximation. This accommodates the com- 
pactness of the superfluid phase of the ordering fields and ac- 
counts properly for the current-current interaction. The Vil- 
lain Hamiltonian for the two-component condensate is given 
by 



-/5/2 Pi(xi-2Tmi ,j) +P2(X2-27rn2,fi)' -Pd (Xl -X2 -27r(ni_ -na,^,)) 



(22) 



r 



where the partition function of the system is given by Z = 
jl"" VeiV92e-'f^"\ and (3 = l/ksT. We have performed 
Monte Carlo calculations on Eq. ( |22] |. using local Metropo- 
lis updating of the fields, 9i{r),92{r) G [0, 27r), while en- 
suring that A6'j(r) G [— 7r,7r). The system sizes considered 
were L x L x L with L = 16, 24, 32, 40, 48, 56 and 64. We 
have chosen /3 = 1 and varied p ~ pi — p2/a. Addition- 
ally, the drag pd is chosen proportional to p, and thus, there 
is technically no difference between this approach and vary- 
ing the temperature for fixed p, pd- During the computations, 
we sample the total energy Hy of the system, and various 
helicity moduli. There are six different helicity moduli we 
keep track of (not all independent). The most general helic- 
ity modulus one can define in this system is applying a twist 
9i ^ 9i + air ■ e^J and 6*2 — * 6*2 + a2r ■ 6^(5. The he- 
licity modulus is then given as the second derivative of the 
free energy with respect to S. For details, see Appendix |A] 
We measure the helicity modulus associated with six differ- 
ent choices of twists, (ai — l,a2 — 0), (fli = 0,02 = 1), 
(ai = 1,02 = ±1) and (ai — l,a2 ~ ±2) i.e. twists in 
01, ^2, 9i ± 02 and 9i ± 26'2, respectively. These are denoted 



T^', T^, and T'^^^^. Here, = T^' ± 2X5^2 



T2 and 



1,±2 



AT 2. A finite helicity modulus is a 
signal of a finite superfluid density of the associated quantity, 
a finite represents the possibility of having co-(counter- 
)superflow of the two components. Likewise, the vanishing of 
the helicity moduli T(^^ signals a thermally driven sponta- 
neous proliferation (blowout) of vortex loops originating with 
multiples of 27r-windings in the phases ai9i + 026*2- We have 
considered these quantities for equal as well as for different 
bare phase stiffnesses pi and p2, and have in all cases varied 
the drag coefficient pd from up to the maximum allowed 
value compatible with the stability of the two-component su- 
perfluid ground state. The location of the phase transitions are 
read off from the peak in the heat capacity 



We first discuss the case pi = p2, for which results for 
the phase diagram and helicity moduli are shown in Fig. |4] 
The dotted lines represent the predictions based on our ana- 
lytical arguments from the previous sections. At pd = 0, the 
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FIG. 4: (Color online) Phase diagram and a set of helicity moduli 
for the model Eq. i22l with equal bare stiffnesses . The shaded 
region illustrates the forbidden parameter regime pd > p/2. The 
helicity moduli are Ti,T2, and T_. The leftmost helicity moduli 
are measured for a drag pd = 0.30p, while the rightmost for pd ~ 
0.39/9. 

system features a doubly degenerate phase transition from a 
2-component superfluid to a 2-component normal fluid at the 
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critical couplings pci = Pc2 = 0.33. These phase transitions 
are in the 3DXy-universality class. When drag is introduced, 
it initially has the effect of reducing the stiffnesses of the in- 
dividual phases Oi and 62, thus moving the doubly degener- 
ate phase transitions to higher couplings (pic, P2c)- At large 
enough drag these phase transitions split, and the intermedi- 
ate phase with ordering only in the phase sum emerges (the 
"paired superfluid phase" in terms of Ref. 3). We observe that 
our computations show that the analytic arguments advanced 
in previous sections describe quite accurately the phase dia- 
gram. 

The line of transition from J7(l)x[/(l)toa fully symmetric 
phase changes its slope indicating that composite vortices for 
sufficiently large drag initiate the transition into the symmet- 
ric state (the preemptive vortex-loop proliferation scenario). 
Importantly, near the bending point the actual transition line 
is situated to the right of the dotted lines, which originate with 
the above bare-stiffness arguments when sub-leading type of 
topological defects are not taken into account. Therefore, 
these estimates naturally underestimate the stiffness at the ac- 
tual position of a preemptive transition. However, even in this 
region, the deviation is not significant. 

The transition line from the state with ordering only in the 
phase sum to a fully symmetric state precisely coincides with 
the analytic estimates and is independent of pd, in the equal 
stiffnesses case, away from the splitting point. The splitting 
point takes place at significantly higher coupling constants in 
the phase diagram than what the naive energy-scale based ar- 
gument gives, and is also in good agreement with the splitting 
point of the preemptive loops proliferation scenario discussed 
in Sections HVCl and IVAl 

The corresponding results for the various helicity moduli 
are also shown in Fig. |4] In the lower right panel the helicity 
modulus T_ for the composite vortex mode (1,-1) vanishes 
first as we approach lower couplings (or equivalently, higher 
temperatures) from the completely ordered side. The resulting 
state is only partially ordered. The individual stiffnesses Ti 
and T2 vanish simultaneously at some lower coupling (higher 
temperature), rendering the system a normal fluid. The in- 
teresting part of the phase diagram is just below the splitting 
point, where we have a region in which the phase transition is 
first-order. 



We find strong indications, shown in Fig. |5] that the transi- 
tion from the [/(I) x J7(l) state to the fully symmetric state in 
the region where vortex-matter based argument suggest pre- 
emptive scenario is indeed a first order transition. This is 
also in agreement with previous computations of the J-current 
model^. 

We proceed to discuss the case of slightly unequal bare 
stiffnesses, i.e. (1,— n) vortices with n > 1 are unimpor- 
tant. In our computations, we have used p2 = l.lpi, see 
Fig. |6] At = the system features two independent phase 
transitions in the BOXF-universality class at pd ~ 0.33 and 
Pc2 ~ 0.30. When drag is introduced, it initially has the effect 
of driving the transitions to higher values of p (lower values 
of T). For moderate values of drag, these two transition close 
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FIG. 5: (Color online) The energy histograms for {p, pd) ~ 
(0.60, 0.20) , with a — 1, i.e. in the preemptive region. A clear 
double peak structure is seen to develop, an indication of a first order 
transition. The areas under the histograms are normalized to 1. 



in on each other, before they merge into one transition from a 
[/(I) X U{1) state into the symmetric state. In terms of vortex 
matter, this is the preemptive region of the phase diagram. For 
even larger drag this line splits, and the intermediate phase 
with ordering associated with the phase sum emerges. The 
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FIG. 6: (Color online) Phase diagram and a set of helicity moduli 
for the model Eq. l l22t . for a — 1.1. The shaded region illustrates 
the forbidden parameter regime pd > p\pil{p\ p2). The helicity 
moduli are Ti,T2, and T_. The left most helicity moduli are mea- 
sured for a drag pd = 0.25/5, while the rightmost are pd = 0.39p. 

dotted lines in Fig. |6]are predictions described in Section IVl 



and these agree well with our computations. Specifically, we 
observe that when the helicity modulus T_ is renormalized 
to zero, the individual stiffnesses become equal, as expected 
from our separation of variables arguments Eq. (fl4] l. More- 
over, in Fig. |7] we show the corresponding energy histograms 
computed on the phase-transition line between points (b) and 
(c) in Fig. |6] namely at {p,pd) ~ (0.60,0.22). This puts us 
in a part of the phase diagram where we would expect, based 
on our vortex-matter arguments, to be able to see the preemp- 
tive scenario explained above played out. Indeed, the phase 
transition is clearly seen to be of first order also in this case, 
thus confirming that the preemptive vortex-loop proliferation 
scenario is also realized for unequal bare phase stiffnesses pi 
and p2- 
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FIG. 7: (Color online) The energy histograms for (p, pd) ~ 
(0.60, 0.22), with a — 1.1, i.e. in the preemptive region. A clear 
double peak structure is seen to develop, an indication of a first order 
phase transition. The areas under the histograms are normalized to 
1. 



We now discuss the case of significantly different bare stiff- 
nesses, i.e. when (1, — n)-vortices with n > 1 are important. 
In our computation we have used p2 = O.SSpi, which from 
Table H] indicate that we should observe a state with order in 
the weighted phase sum, with n = 2. As in the case of slightly 
unequal stiffnesses the system features two independent tran- 
sition in the 3DXy-universality class, in our computation the 
transitions at pd ^ occurs at pd w 0.33 and pc2 ~ 0.605. 
At small drag values the transitions stay independent and are 
shifted to higher values of p. For moderate drag values the 
region with partial order (order in 6i) becomes smaller, be- 
fore disappearing at some higher drag value. The system then 
enters into the preemptive vortex-loop proliferation region, 
where the system features a transition from a U{1) x U{1)- 
state to the fully symmetric state. In this region we find strong 
indications of first order transitions, shown in Fig. |9] 
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FIG. 8: (Color online) Phase diagram and a set of helicity moduli 
for the model Eq. i22i . a = 0.55. The shaded region illustrates 
the forbidden parameter regime pd > pip2/{pi + P2). The helicity 
moduli are Ti,T2, and Ti,-2. Here, the latter correspond to T_, 
with the difference that the -phase is twisted twice as much as 61. 
The leftmost helicity moduli are measured for a drag pd = 0.32p, 
while the rightmost are pd ~ 0.336p. 
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FIG. 9: (Color online) The energy histograms for (p, pd) ~ 
(1.77, 0.58), with a — 0.55, i.e. in the preemptive region. A clear 
double peak structure is seen to develop, an indication of a first order 
transition. The areas under the histograms are normalized to 1 . 



For even larger drag values this single lines splits into two 
lines, and a partially ordered state appears. The partially or- 
dered state which appears is a state were (1, — 2)-vortices have 
proliferated while individual vortices stay confined (in gen- 
eral (1, — n)-vortices can proliferate). In the lower rightmost 
panel in Fig. [8]we observe that Ti _2 drops to zero while Ti 
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and T2 remain finite, from this we conclude that we observe 
the partially ordered state with order in the weighted phase 
sum. The MC calculation shows that in this quite generic case 
that the analytical arguments from the first part of the paper 
are remarkably quantitatively accurate even near line-splitting 
points. Observe further that when Ti _2 have renormalized 
to zero we have the relation Ti = 4T2, as expected from the 
discussion nearEq. (l20b (in general we expect Ti = n^T2). 



VIII. IMPLICATION FOR THE LIQUID METALLIC 
HYDROGEN PROBLEM 

The approach developed above is also useful to obtain 
insight into the role of non-dissipative drag when it is in- 
cluded in the problem of multicomponent electrically charged 
condensates^iiS. For example, in the problem of projected 
quantum fluid states of hydrogen, we deal with two electri- 
cally charged fields corresponding to electronic and protonic 
condensates. The fields will then be coupled by an electro- 
magnetic gauge field in addition to the now familiar drag cou- 
pling 

F = \ I dr{p(V6'i -eAf + ap{Ve2 + eAf 

-Pd (V^^l - V02 - 2eA)2 + (V X A)' }, 

(23) 

were e is the charge and A is the gauge field. 

In this case, the physically relevant separation of variables 
corresponds to extraction of the phase sum. This follows, 
since in this situation it is the phase sum which is not cou- 
pled to the gauge field. This allows us to draw conclusions 
about superfluid and superconducting states of the system^. 
Following the same line of reasoning, an assessment of the 
role of non-dissipative drag is made by an extraction of the 
phase sum to distinguish the drastically different charged and 
neutral modes of the system. The model then becomes 



1 



■ [{p - 2pd) V01 - {ap ~ 2pd) V62 



(1 + a)p - Apd 

-e{{l + a)p-Apd}AY + {V X Afy (24) 



By virtue of featuring one composite charged mode and one 
composite neutral mode, this model has the same structure 
as the model with zero drag^ '". However, now the stiff- 
nesses of neutral and charged modes acquire dependence to 
the drag coefficient pd- Therefore, the conclusions of Ref. |^ 
should be rather robust against finite-drag perturbations. The 
inter-component drag term is a quite different perturbation 
to the system compared to inter-component Josephson cou- 
pling. The latter is prohibited in hydrogen, but is allowed in 
multicomponent electronic condensates. Josephson coupling 
amounts to an explicit symmetry breakdown, and in terms of 
long length scale physics it represents a singular perturbation 



compared to the case where it is absent. An inter-component 
drag term has two gradients in it, since it is a current-current 
interaction. Consequently, it has a naive scaling dimension 
which is reduced by 2 compared to the Josephson coupling, 
and in contrast to the Josephson couplling it does not repre- 
sent a singular perturbation. Quite the contrary, as we have 
seen, a critical value of the strength of the inter-component 
drag term is required for it to have an appreciable effect on 
the physics of the system. 



IX. SUMMARY AND CONCLUSIONS 

In this paper, we have studied the problem of the influence 
of non-dissipative inter-component drag on the phase diagram 
and phase transitions in a two-component Bose-Einstein con- 
densate. The non-dissipative drag is a quite generic feature 
present in interacting multicomponent systems in the contin- 
uum as well as on a lattic o'i^i^i'^ . Recently, the topology of 
the phase diagram and orders of the phase transitions were in- 
tensively studied in the J-current model with U{1) x U{1) 
symmetry by means of worm-algorithm based Monte Carlo 
simulations^'^-^-^, revealing novel features such as conversions 
of the phase transitions from continuous to first order as a 
function of drag strength. 

We have developed an approach in terms of topological de- 
fects for understanding these phase transitions and get new 
insight into physics of the various states of two-component 
Bose-Einstein condensates. We have carried out an investi- 
gation of the phase diagram based on analytical vortex-matter 
arguments, and suggested a novel scenario of vortex-matter 
behavior, namely a "preemptive vortex-loop proliferation". 
Such a scenario may well be generic to systems where sym- 
metry is restored through proliferation of distinct topological 
defects in the form of vortex loops that have been excited out 
of the individually conserved condensates. We have found 
support for these scenarios in large-scale Monte Carlo calcu- 
lations. These computations have been carried out using a 
representation of the system in terms of the phase of the com- 
plex ordering field of each of the components. The approach 
allows us to investigate directly the physics of topological de- 
fects in this system. Importantly, the phase representation also 
allows us to study the system under rotation. This can provide 
a bridge for studying these states of matter experimentally via 
rotational response. Work on this problem is in progresses. 



Acknowledgments 

The authors acknowledge useful discussions with J. Hove, 
D. A. Huse, A. Kuklov, E.J. Muefler B. Svistunov, and M. 
Wallin. This work was supported by the Norwegian Research 
Council Grant Nos. 1585187/431, 158547/431 (NANOMAT), 
and Grant No. 167498/V30 (STORFORSK). The authors 
thank the Center for Advanced Study at the Norwegian 
Academy of Science and Letters, where part of this work was 
done. 



11 



APPENDIX A: SUPERFLUID DENSITY MATRIX IN A 
2-COMPONENT SYSTEM 

In general, the helicity modulus defines the superfluid den- 
sity of a system. For the Andreev-Bashkin probleroii, the 
superfluid density is a matrix quantity given by the second 
derivative of the free energy of the system with respect to an 
infinitesimal twist in the phase, i.e. d{r) 9{r) — S r. The 

hehcity modulus, Y, is then given as — -p- ^gf' 

Since F[S] = —(i^^ \n.Z[S\, where j3 is inverse temperature 
andZ[5] = J 2?re~''^[''l is the partition function, the helicity 
modulus can further be written as 



.5=0 



This is a general expression for the helicity modulus, indepen- 
dent of the form of the Hamiltonian. We now specify the form 
of the Hamiltonian to that of a two-component Villain-model 
i.e. = X;^_^K(A,.6li(r), Ai.6'2(r)) where the potential 
is given in Eq.(l22]l. We now apply an arbitrary twist in 

the phases, (^^j^j) ^ (^^[^j )- (;^^ ) <5 • r, were ai, are 

two real numbers and expressions on both side of the arrow 
satisfies periodic boundary conditions. The Hamiltonian then 
takes the form [8] ^ J2rM K ( A.^i {r)-aiS,, A,02 (r) - 
02 The first and second derivatives of the Hamiltonian are 
then given by. 




dH[5] 



d''H[S] 



5=0 



E 

r 



-ai 



2 d'V, 



0,2- 



dAJ, 'dAJ2 



2aia2 



6=0 r ^ -"M-l ^^^..i^^p. 

The helicity modulus associated with this choice of twist in the phase, is given by 



9A^^2 



2aia2 



L3 




dA^ei 




dA^e^dA^d2 



dA.e2 



^dA^ei 



( dH, 



\dA^e2 \dA^92 



(A2) 
(A3) 

(A4) 
(A5) 
(A6) 



We observe that a general twist in the phases can be expressed through three independent quantities, the superfluid density of 
the two single components eq. ( IA4b and iA6h denoted T^* and Tt^ respectively, and a novel inter-component quantity Eq. 
iA5\ denoted T^2- We interpret T12 as a renormalized drag coefficient. A general helicity modulus may then be written in the 
compact form 



ai ,02 



afTi + 2aia2T]'2 + 0,2^2. 

I 



(A7) 



APPENDIX B: NEGATIVE DRAG COEFFICIENT 

The subject of the sign of the drag coefficient, pd is a subtle 
one and depends on the physical realization of the model. In 
the case of a realization of the model on an optical lattice the 
sign of the drag coefficient can straightforwardly be made 
negative^. 

In the case of a negative drag coefficient the analysis in Sec- 
tionlVlwill hold, with the role of (1, 1)- and (1, — l)-vortices 
interchanged. However in the separation of variables we need 
to extract a phase difference to estimate the stiffness which 



would remain in the system when (1. 1) vortices proliferate. 
Then, the proper separation of variables for analyzing the 
model is given by 



F = dr 



l + a 
P 

l + a 



" p-pd][V{0i-02y' 



[V6I1 + aV6l2] 



(Bl) 
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FIG. 10: (Color online) Same as Fig. (2] but now also including a 
negative drag coefficient pd- In region I) the system features a bro- 
ken (/(I) X (7(1) symmetry, in regions II), IV) and V) the system 
features a broken (7(1) symmetry associated with the phase sum, one 
individual phase and the phase difference respectively. In region III) 
the system is in the fully symmetric state. (See also Fig. 
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